\chapter{Reconstruction of Parallel MRI} \label{chap:4}
\section{Introduction}
The concept of parallel MRI was first suggested by Carlson \cite{Carlson1987}, but his work was unknown until 2004 \cite{Griswold2007}. The introduction of simultaneous acquisition of spatial harmonics (SMASH) \cite{Sodickson2001} started the age of parallel MRI.

The conventional Fourier MRI is modeled by
\begin{align*}
    y(t) = \int M(\mathbf{x}) e^{\imath \mathbf{k}(t) \cdot \mathbf{x}} d \mathbf{x}
\end{align*}

There is more than one RF coil in a parallel MRI system. Each coil has a different receiving property that is called localized sensitivity. Let $C_{k}(\mathbf{x})$ be the localized sensitivity of the $k$th coil, $k=1,\cdots,K$. The received signal at the $k$th coil can be modeled by
\begin{align}
    y_k(t) = \int C_k (\mathbf{x}) M(\mathbf{x}) e^{\imath \mathbf{k}(t) \cdot \mathbf{x}} d \mathbf{x}
    + \varepsilon _k (t) \label{eq:pMRI}
\end{align}
where $\varepsilon _k (t)$ is the additive white Gaussian noise with zero mean and $\sigma _{\varepsilon}^2$ variance.

\cite{DIETRICH2007a} summarized the advantages of parallel imaging:
\begin{itemize}
    \item Faster imaging.
    \item Higher spatial resolution.
    \item Improved image quality of single-shot or turbo-spin-echo or echo-planar by shortening the echo train which is severely affected by signal decay and field inhomogeneities.
    \item Parallel imaging ideally complements the increased SNR and compensates the growing specific absorption rate and increasing geometric distortions for high and ultra-high field MRI.
\end{itemize}

The SNR (signal-to-noise ratio) can be increased with array coils in the same imaging time. The array coils can also be used for partially parallel acquisition (PPA) \cite{Griswold2000}. In this chapter, we extend the parallel MRI model \eqref{eq:pMRI} with the concept of SS-PARSE.

In current reconstruction algorithms, the coil sensitivities $C_k(\mathbf{x})$ are found by putting an known object in a parallel MRI system and measuring the spatial sensitivities directly. This method is convenient, but it is inaccurate because different objects may have different coil sensitivities due to the fact that each object has its own properties that interact with the coil responses. We can bring the concept of SS-PARSE into parallel MRI. We call this extension parallel SS-PARSE (PSSPARSE). The received signal in the $k$th coil is modeled by:
\begin{align}
    y_k(t) = \int C_k (\mathbf{x}) M_0(\mathbf{x}) e^{\left[ R_2^* (\mathbf{x}) + \imath \omega (\mathbf{x}) \right]t} e^{\imath \mathbf{k}(t) \cdot \mathbf{x}} d \mathbf{x}
    + \varepsilon _k (t) \label{eq:pPARSE}
\end{align}
The goal of PSSPARSE is to reconstruct $C_{k}(\mathbf{x})$, $M_{0}(\mathbf{x})$, $R_{2}^{*}(\mathbf{x})$ and $\omega(\mathbf{x})$ from the observed signals $y_k (t)$ for all $k=1,\cdots,K$.

\section{Reconstruction}
\subsection{Extension of SS-PARSE} \label{sec:4-2-1}
By discretizing \eqref{eq:pPARSE} on the spatial $(x,y)$ grid indexed by $i$, we have
\begin{align}
    y_k (n) = \underbrace{\sum _i C_{k,i} M_{0i} e^ {nW_i} e^{\imath \mathbf{k}_n \cdot \mathbf{x}_i}}
    _{\hat{y}_k}+ \varepsilon _k (n) \label{eq:Dis-pPARSE}
\end{align}
where $n=1,\cdots,N$, $W_i = -\left[ R_{2i}^{*} + \imath \omega_{i} \right] \Delta t$, $\mathbf{k}_{n} = \mathbf{k}(n\Delta t)$, and $\Delta t$ is the sampling interval.

We use the method discussed in Chapter \ref{chap:2} to simultaneously solve $C_k$, $M_0$ and $W$. We use iterative conjugate-gradient algorithm to minimize the cost function:
\begin{align}
    J(\mathbf{z}) = \sum_{k=1}^{K} \left\| \mathbf{y}_k - \hat{\mathbf{y}}_k (\mathbf{z}) \right\| ^{2} \label{eq:pCost}
\end{align}
with respect to $\mathbf{z}$, where $\mathbf{z}=\{\mathbf{C}_k,\mathbf{M}_0, \mathbf{W}, k=1,\cdots,K\}$.

The gradient $\nabla _{\mathbf{z}}J$ computations include three parts:
\begin{align}
    \frac{\partial J}{\partial \mathbf{M}_0} &= \sum_{k=1}^{K} \mathbf{C}_k
    \sum _{n=1}^{N} f_{k}(n) e^{n\mathbf{W}} e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pGrad-M} \\
    \
    \frac{\partial J}{\partial \mathbf{W}} &= \mathbf{M}_0 \sum_{k=1}^{K} \mathbf{C}_k
    \sum _{n=1}^{N} nf_{k}(n) e^{n\mathbf{W}} e^{\imath \mathbf{k}_{n} \cdot \mathbf{x}} \label{eq:pGrad-W} \\
    \
    \frac{\partial J}{\partial \mathbf{C}_k} &= \mathbf{M}_0 \sum _{n=1}^{N} f_{k}(n) e^{n\mathbf{W}} e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pGrad-S}
\end{align}
where $f_{k}(n) =[ \hat{y}_{l} (n) -y_{k}(n)]^{*}$.

We apply the polynomial approximation of the exponential time function in Chapter \ref{chap:2} to compute the cost function and gradient for parallel MRI.
\begin{align}
    y_k(n) & \approx e^{jn\omega _0} \sum _{l=0}^{L-1} n^{l} \sum _{i=1}^{K} C_{k,i}M_{0i} Z_{i} (l) e^{\imath \mathbf{k}_n \cdot \mathbf{x}(i)} \label{eq:Fast-Dis-pPARSE} \\
    \
    \frac{\partial J}{\partial \mathbf{M}_0} &\approx \sum_{k=1}^{K} \mathbf{C}_k
    \sum_{l=0}^{L-1} \mathbf{Z} (l) \sum _{n=1}^{N} f_{k}(n) e^{ \imath n \omega_0}n^{l}
    e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pFast-Grad-M} \\
    \
    \frac{\partial J}{\partial \mathbf{W}} &\approx \mathbf{M}_0 \sum_{k=1}^{K} \mathbf{C}_k
    \sum_{l=1}^{L} \mathbf{Z}(l) \sum _{n=1}^{N} f_k(n) e^{\imath n\omega_0}n^{l}
    e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pFast-Grad-W} \\
    \
    \frac{\partial J}{\partial \mathbf{C}_k} &\approx \sum_{l=0}^{L-1} \mathbf{Z} (l)
    \sum _{n=1}^{N} f_{k}(n) e^{ \imath n \omega_0}n^{l}
    e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pFast-Grad-C}
\end{align}
$\mathbf{Z}$ and $\omega _{0}$ are defined in chapter \ref{chap:2}.

The reconstruction algorithm based on cubic interpolation can also be applied to the reconstruction in parallel MRI.

\subsection{Initialization}
The first step in the reconstruction of PSSPARSE is to find the initial values of $C_k(\mathbf{x})$, $M_{0}(\mathbf{x})$, $R_{2}^{*}(\mathbf{x})$ and $\omega(\mathbf{x})$ for the conjugate-gradients algorithms described above. We use a different mathematical model to find the initial values. We assume that there are $K$ different $M_0$ maps --- $M_1, \cdots,M_K$ --- the magnitude. Based on this assumption, $y_k$ is modeled by:
\begin{align}
    y_k (n) = \underbrace{\sum _i M_{k,i} e^ {nW_i} e^{\imath \mathbf{k}_n \cdot \mathbf{x}_i}}
    _{\hat{y}_k}+ \varepsilon _k (n) \label{eq:Multi-M0-pPARSE}
\end{align}
In the cost function \eqref{eq:pCost}, the unknown $\mathbf{z}$ is defined as $\mathbf{z}=\{\mathbf{M}_k, \mathbf{W}, k=1,\cdots,K\}$.

The gradients are computed by:
\begin{align}
    \frac{\partial J}{\partial \mathbf{M}_k} &= \sum _{n=1}^{N} f_{k}(n) e^{n\mathbf{W}}
    e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:Multi-M0-Grad-M} \\
    \
    \frac{\partial J}{\partial \mathbf{W}} &= \sum_{k=1}^{K} \mathbf{M}_k
    \sum _{n=1}^{N} nf_{k}(n) e^{n\mathbf{W}} e^{\imath \mathbf{k}_{n} \cdot \mathbf{x}} \label{eq:Multi-M0-Grad-W}
\end{align}

With polynomial approximation, the estimated signals and and gradients are evaluated by:
\begin{align}
    y_k(n) & \approx e^{jn\omega _{0}} \sum _{l=0}^{L-1} n^{l} \sum _{i=1}^{K} M_{k,i}
    Z_{i} (l) e^{\imath \mathbf{k}_{n} \cdot \mathbf{x}(i)} \label{eq:Fast-Multi-M0-PARSE} \\
    \
    \frac{\partial J}{\partial \mathbf{M}_k} &\approx \sum_{l=0}^{L-1}
    \mathbf{Z} (l) \sum _{n=1}^{N} f_{k}(n) e^{ \imath n \omega_0}n^{l}
    e^{\imath \mathbf{k}_{n} \cdot \mathbf{x}} \label{eq:Fast-Multi-M0-Grad-M} \\
    \
    \frac{\partial J}{\partial \mathbf{W}} &\approx \sum_{k=1}^{K} \mathbf{M}_k
    \sum_{l=1}^{L} \mathbf{Z}(l) \sum _{n=1}^{N} f_k(n) e^{\imath n\omega_{0}}n^{l}
    e^{\imath \mathbf{k}_n \cdot \mathbf{x}} \label{eq:pFast-Multi-M0-Grad-W}
\end{align}

After finding $\mathbf{M}_k$, we use the root mean square of $\mathbf{M}_k$ as the initial values of $\mathbf{M}_0$:
\begin{align}
    \mathbf{M} _0 = \sqrt{\frac{\sum _{k=1}^{K} \left| \mathbf{M} _k \right| ^2} {K}} \label{eq:M0-Est}
\end{align}
and the coil sensitivities are computed by:
\begin{align}
    \mathbf{C} _k = \frac{\mathbf{M}_k}{\mathbf{M}_0},~~k=1,\cdots,K \label{eq:Coil-Est}
\end{align}

We also apply the interpolation method to compute the initial conditions. Because of the use of interpolation, the algorithm stated in \ref{sec:4-2-1} is initialized with the proper interpolation coefficients that are associated with $\mathbf{M}_0$ and $\mathbf{C}_k$.

\subsection{Regularization}
We use regularization to improve the reconstruction performance. The regularization operation is defined as:
\begin{align}
    R(\mathbf{x}) = \mathbf{H}*\mathbf{x}*\mathbf{H}^T \label{eq:Reg-Op}
\end{align}
where $*$ denotes convolution operation. The regularization kernel $\mathbf{H}$ is:
\begin{align}
    \mathbf{H} = \begin{bmatrix}
    0 & 0 & 0 \\ -1 & 2 & -1 \\ 0 & 0 & 0
    \end{bmatrix} \label{eq:reg-kernel}
\end{align}

With the regularization terms, the cost function \eqref{eq:pCost} is changed to:
\begin{align}
    J(\mathbf{z}) = & \sum_{k=1}^{K} \left\| \mathbf{y}_k - \hat{\mathbf{y}}_k
    (\mathbf{z}) \right\| ^{2} \notag \\
    &+ \alpha \left\| R(\mathbf{M}_0) \right\| ^2
    + \beta \sum _{k=1}^{K} \left\| R(\mathbf{C}_k) \right\| ^2
    + \gamma _R \left\| R(\mathbf{R}_2^*) \right\| ^2
    + \gamma _I \left\| R(\bbomega) \right\| ^2 \label{eq:Reg-pCost}
\end{align}
where $\alpha$, $\beta$, $\gamma _R$ and $\gamma _I$ are all nonnegative real numbers.

The corresponding gradients $\nabla J$ includes terms of $\nabla R$:
\begin{align}
    \frac{1}{2} \nabla \left\| R (\mathbf{x}) \right\| ^2 = \mathbf{A} * \mathbf{x} \label{eq:Grad-Reg}
\end{align}
where
\begin{align}
    \mathbf{A} = \begin{bmatrix}
    0 & 0 & 1 & 0 & 0 \\
    0 & 0 & -4 & 0 & 0 \\
    1 & -4 & 12 & -4 & 1 \\
    0 & 0 & -4 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0
    \end{bmatrix} \label{eq:Grad-Reg-Matrix}
\end{align}

\section{Simulation}
We synthesized the simulation data from the images in Figure \ref{fig:fig-Std-Objs} and \ref{fig:fig-Std-Coils}. Eq. \eqref{eq:Dis-pPARSE} was used to to generate the synthesized data. The image resolution for the synthesis was $1024 \time 1024$.

It is assumed that field of the view is 12.8cm. Both $x$-axis and $y$-axis are defined on $[-6.4,6.4]$. Four receiving coils are placed at the four corners of the FOV. The coil sensitivities are described by \eqref{eq:Sim-Coil-Sen}.
\begin{align}
    C(x,y) &= \left[ (x \pm x_0)^2 + (y \pm y_0)^2 \right] ^{-\alpha}
    \label{eq:Sim-Coil-Sen}
\end{align}
In the simulation, $x_0 = y_0 = 6.8$ and $\alpha = 1/4$.

The SNR is defined as:
\begin{align}
    \text{SNR} = \frac{\sum _{k=1}^{K} \left\| \hat{\mathbf{y}}_k \right\| ^2}
    {MN \sigma ^2 _\varepsilon}
\end{align}

\begin{figure*}
    \centering
    \subfloat[$|M_{0}|$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_M0.eps}}
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_R2e.eps}} \\
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_Frmap.eps}} \\
    \caption{The magnitude, decay and field map used to synthesize simulation data. All images are displayed with $256 \times 256$ resolution. $M_0$ is normalized to $1$.}
    \label{fig:fig-Std-Objs}
\end{figure*}
\begin{figure*}
    \centering
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_C1.eps}}
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_C2.eps}} \\
    \subfloat[Coil 3]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_C3.eps}}
    \subfloat[Coil 4]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Std_C4.eps}} \\
    \caption{The coil sensitivity ,maps used to synthesize simulation data. All images are displayed with $256 \times 256$ resolution. The unconstrained area in which $M_0$ is zero are displayed with zero. All maps are normalized to $1$.}
    \label{fig:fig-Std-Coils}
\end{figure*}

\subsection{Interpolation of Coil Sensitivity}
In the all of the simulations, $128 \times 128$ is used as image resolution. We use $2 \times$ interpolation for the computation of the initial conditions. That is, the coefficients of $\mathbf{M}_k$ and $\mathbf{W}$ are $64 \times 64$. In the joint reconstruction of  $\mathbf{M}_0$, $\mathbf{W}$ and $\mathbf{C}_k$, $64 \times 64$ resolution is also the choice of the coefficients of $\mathbf{M}_0$ and $\mathbf{W}$. Since the $\mathbf{C}_k$'s are very spatially smooth, we use a lower resolution for the coefficients of $\mathbf{C}_k$'s --- a larger interpolation factor. We used the cubic convolution interpolation for all of the estimated parameters. For $I \times$ interpolation, the sampling interval is $h=1/I$. The 1-D interpolation vector of $I \times$ interpolation is:
\begin{align}
    u( \pm n h), ~~~n= 0,\pm h, \cdots, \pm (3/h-1)h
\end{align}
where $u(x)$ is defined by \eqref{eq:CubicIntrpKrl}.

We experimentally compared the reconstruction accuracies of different interpolation factors of the coil sensitivity in Table \ref{tab:tab-Coil-Intrp-Factor}. In this comparison, we used signals with 30 dB SNR.
\begin{table}[h!]
\centering
\caption{NRMSE (\%) of Different Coil Sensitivity Interpolations}
\begin{tabular}{cccccccc}
    Interpolation & $C_1$ & $C_2$ & $C_3$ & $C_4$ & $M_{0}$ & $R_{2}^{*}$ & $\omega$ \\
    \midrule
    $2 \times$ & 29.2 & 36.8 & 36.6 & 30.3 & 29.8 & 14.4 & 17.9 \\
    \midrule
    $4 \times$ & 8.7 & 8.7 & 7.5 & 7.2 & 28.3 & 14.4 & 17.9 \\
    \midrule
    $8 \times$ & 6.9 & 8.3 & 7.5 & 6.7 & 29.5 & 14.4 & 17.9 \\
    \midrule
    $16 \times$ & 10.7 & 10.4 & 6.9 & 7.2 & 30.2 & 14.4 & 17.9 \\
    \bottomrule
\end{tabular}
\label{tab:tab-Coil-Intrp-Factor}
\end{table}

The reconstructed $C_1$ from the different interpolation factors are illustrated in Figure \ref{fig:fig-Coil-Intrps}.
\begin{figure*}
    \centering
    \subfloat[$2 \times$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_02x_C1_30dB.eps}}
    \subfloat[$4 \times$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_04x_C1_30dB.eps}} \\
    \subfloat[$8 \times$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_08x_C1_30dB.eps}}
    \subfloat[$16 \times$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_16x_C1_30dB.eps}} \\
    \caption{The coil sensitivity $C_1$ reconstructed from different interpolation factors. All images are displayed with $128 \times 128$ resolution. The artifacts in the unconstrained area in which $M_0$ is zero are removed. All maps are normalized to $1$. No regularization is used.}
    \label{fig:fig-Coil-Intrps}
\end{figure*}

\begin{figure*}
    \centering
    \subfloat[$|M_{0}|$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_M0.eps}}
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_R2e.eps}} \\
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_Frmap.eps}} \\
    \caption{The magnitude, decay and field map reconstructed from signals with 30dB SNR. All images are displayed with $128 \times 128$ resolution. The artifacts in the unconstrained area in which $M_0$ is zero are removed. $M_0$ is normalized to $1$. No regularization is used.}
    \label{fig:fig-30dB-Objs-NoReg}
\end{figure*}
\begin{figure*}
    \centering
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_C1.eps}}
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_C2.eps}} \\
    \subfloat[Coil 3]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_C3.eps}}
    \subfloat[Coil 4]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_203_30dB_C4.eps}} \\
    \caption{The coil sensitivity maps reconstructed from signals with 30dB SNR. All images are displayed with $128 \times 128$ resolution. The artifacts in the unconstrained area in which $M_0$ is zero are removed. All maps are normalized to $1$. No regularization is used.}
    \label{fig:fig-30dB-Coils-NoReg}
\end{figure*}

\subsection{Regularized Reconstruction}
We applied the regularization to this reconstruction. Because the field map $\omega$ is not spatially smooth enough, the experiments show that the regularization hardly improve the reconstruction accuracy of $\omega$. We empirically select the regularization parameters $\alpha$, $\beta$ and $\gamma _R$. Table \ref{tab:tab-Regularization} compares the reconstruction performance for different combinations of regularization coefficients. The interpolation factor for the coil sensitivity used here is $8$. The interpolation factors of $M_0$, $R_2^*$ and $\omega$ are all 2. A set of coil sensitivities and images from regularized reconstruction are displayed in Figure \ref{fig:fig-30dB-Objs-Reg} and \ref{fig:fig-30dB-Coils-Reg}. Table \ref{tab:tab-Regularization} shows that the selection of the regularization parameter for one set of variables has little impact on the reconstruction accuracy of the other sets of variables. For example, the selection of $\alpha$, the regularization parameter of $M_0$, is almost unrelated to the reconstruction accuracy of $R_{2}^{*}$ and $\omega$.
\begin{table}[h!]
\centering
\caption{NRMSE (\%) of Different Regularization Coefficients}
\begin{tabular}{cccccccccc}
    $\alpha$ & $\gamma _R$ & $\beta$ & $C_1$ & $C_2$ & $C_3$ & $C_4$ & $M_{0}$ & $R_{2}^{*}$ & $\omega$ \\
    \midrule
    0 & 0 & 0 &
    7.0 & 8.3 & 7.5 & 6.7 & 29.5 & 14.4 & 17.9 \\
    \midrule
    0 & 0 & $1.09 \times 10^6$ &
    7.7 & 7.9 & 6.9 & 7.6 & 29.0 & 14.4 & 17.9 \\
    \midrule
    0 & $1.75 \times 10 ^8$ & 0 &
    8.6 & 6.8 & 6.7 & 9.3 & 28.7 & 4.6 & 17.9 \\
    \midrule
    0 & $1.75 \times 10 ^8$ & $1.09 \times 10^6$ &
    10.6 & 8.9 & 8.2 & 8.5 & 28.8 & 5.5 & 17.9 \\
    \midrule
    $5.5 \times  10^7$ & 0 & 0 &
    10.8 & 8.5 & 6.4 & 9.6 & 13.7 & 14.4 & 17.9 \\
    \midrule
    $5.5 \times  10^7$ & 0 & $1.09 \times 10^6$ &
    9.4 & 8.4 & 6.4 & 9.8 & 13.6 & 14.4 & 17.9 \\
    \midrule
    $5.5 \times  10^7$ & $1.75 \times 10 ^8$ & 0 &
    11.5 & 8.9 & 6.6 & 9.5 & 13.7 & 4.0 & 17.9 \\
    \midrule
    $5.5 \times  10^7$ & $1.75 \times 10 ^8$ & $1.09 \times 10^6$ &
    9.5 & 8.4 & 6.5 & 9.5 & 13.6 & 3.9 & 17.9 \\
    \bottomrule
\end{tabular}
\label{tab:tab-Regularization}
\end{table}

\begin{figure*}
    \centering
    \subfloat[$|M_{0}|$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_M0_30dB_Reg.eps}}
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_R2e_30dB_Reg.eps}} \\
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Frmap_30dB_Reg.eps}} \\
    \caption{The magnitude, decay and field map reconstructed from signals with 30dB SNR. The regularization parameters $\alpha = 5.5 \times 10 ^7$, $\gamma _R = 1.75 \times 10^8$ and $\gamma _I = 0$. All images are displayed with $128 \times 128$ resolution. The artifacts in the unconstrained area in which $M_0$ is zero are removed. $M_0$ is normalized to $1$.}
    \label{fig:fig-30dB-Objs-Reg}
\end{figure*}
\begin{figure*}
    \centering
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_C1_30dB_Reg.eps}}
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_C2_30dB_Reg.eps}} \\
    \subfloat[Coil 3]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_C3_30dB_Reg.eps}}
    \subfloat[Coil 4]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_C4_30dB_Reg.eps}} \\
    \caption{The coil sensitivity maps reconstructed from signals with 30dB SNR. The regularization parameters $\alpha = 1.09 \times 10^6$. All images are displayed with $128 \times 128$ resolution. The artifacts in the unconstrained area in which $M_0$ is zero are removed. All maps are normalized to $1$.}
    \label{fig:fig-30dB-Coils-Reg}
\end{figure*}

\section{Human Experiment}
We applied the reconstruction method to the experiment of a human brain. A Siemens Tim Trio 3T MRI system was used in this experiment. The system is located at the Department of Neuroscience of the Brown University. The major parameters of this experiment are listed in Table \ref{tab:tab-Para-Human}. The rosette trajectory used in the experiment is plotted in Figure. \ref{fig:fig-pHuman-rosette}.
\begin{table}[h!]
\centering
\caption{Experiment Parameters of Human Experiment}
\begin{tabular}{lrc}
    Parameter & Value & Unit \\
    \midrule
    sampling interval ($\Delta t$) & 5.0 & $\mu$s \\
    \midrule
    field of the view (FOV) & 22.0 & cm \\
    \midrule
    readout duration & 56.8 & ms \\
    \midrule
    magnetic field & 3.0 & T \\
    \midrule
    trajectory & rosette & - \\
    \bottomrule
\end{tabular}
\label{tab:tab-Para-Human}
\end{table}
\begin{figure*}
    \centering
    \includegraphics[width=4.0in]{C:/TeX/Thesis/Chapter4/Figures/pHumanRosette.eps}\\
    \caption{Rosette trajectory used in human experiment. Only the first half of the trajectory is plotted.}
    \label{fig:fig-pHuman-rosette}
\end{figure*}

In the reconstruction, we used the signals from two coils. Since the simulation experiments show that $4 \times$ and $8 \times$ interpolation for the coil sensitivities give similar results, we applied this two kinds of interpolation. For both cases, $M_0$, $R_2^*$ and $\omega$ were reconstructed with $2 \times$ interpolation. The recovered images are shown in Figure \ref{fig:fig-Brain-Intrp0408-Objs} and \ref{fig:fig-Brain-Intrp0408-Coil}.

We also interpolate coil sensitivity with cubic spline function. The magnitude, decay and field map are still interpolated by cubic convolution function because of its advantage for the exponential time function showed on Chapter 2. The results are in Figure \ref{fig:fig-Brain-Intrp0408-Spline-Objs} and \ref{fig:fig-Brain-Intrp0408-Spline-Coil}.

The reconstructed images show that the proposed algorithm can produce realistic results. The artifacts on the brinks of coil sensitivity maps, $R_2^*$ and $\omega$ are caused by the zero or near zero $M_0$ values at that locations. These artifacts are outside of the field of interest.

$4 \times$ and $8 \times$ coil sensitivity interpolation reconstructed similar $M_0$, $R_2^*$ and $\omega$, but the the coil sensitivity maps from $8 \times$ interpolation are better. We see some artifacts at the left upper and lower corners. The artifacts may be caused by the locations of the receiving coils because the two coils were located at the right upper and lower corners. Two more coils at the left upper and lower corners may correct this problem. Another possible reason for the artifacts is system bias. By tuning reconstruction parameters, the system bias can be removed.
 
\begin{figure*}
    \centering
    \subfloat[$|M_{0}|$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp04_M0.eps}}
    \subfloat[$|M_{0}|$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp08_M0.eps}} \\
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp04_R2e.eps}}
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp08_R2e.eps}} \\
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp04_Frmap.eps}}
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp08_Frmap.eps}} \\
    \caption{The magnitude, decay and field map reconstructed from the experiment of human brain. Cubic convolution interpolation was used to coil sensitivity. All images are displayed with $128 \times 128$ resolution. (a) (c) (e) are from $4 \times$ interpolation for coil sensitivity. (b) (d) (f) are from $8 \times$ interpolation for coil sensitivity. Most of the artifacts outside of the head are removed. $M_0$ is normalized to $1$. No regularization is used.}
    \label{fig:fig-Brain-Intrp0408-Objs}
\end{figure*}

\begin{figure*}
    \centering
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp04_C1.eps}}
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp08_C1.eps}} \\
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp04_C2.eps}}
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_Intrp08_C2.eps}} \\
    \caption{The coil sensitivity maps reconstructed from the experiment of human brain. Cubic convolution interpolation was used to coil sensitivity. All images are displayed with $128 \times 128$ resolution. (a) (c) are from $4 \times$ interpolation. (b) (d) are from $8 \times$ interpolation. Most of the artifacts outside of the head are removed. All maps are normalized to $1$. No regularization is used.}
    \label{fig:fig-Brain-Intrp0408-Coil}
\end{figure*}

\begin{figure*}
    \centering
    \subfloat[$|M_{0}|$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_04x_Spline_M0.eps}}
    \subfloat[$|M_{0}|$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_08x_Spline_M0.eps}} \\
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_04x_Spline_R2e.eps}}
    \subfloat[$R_{2}^{*} (\text{sec}^{-1})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_08x_Spline_R2e.eps}} \\
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_04x_Spline_Frmap.eps}}
    \subfloat[$\omega (\text{Hz})$]{\includegraphics[width=2.5in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_08x_Spline_Frmap.eps}} \\
    \caption{The magnitude, decay and field map reconstructed from the experiment of human brain. Cubic spline interpolation was used to coil sensitivity. All images are displayed with $128 \times 128$ resolution. (a) (c) (e) are from $4 \times$ interpolation for coil sensitivity. (b) (d) (f) are from $8 \times$ interpolation for coil sensitivity. Most of the artifacts outside of the head are removed. $M_0$ is normalized to $1$. No regularization is used.}
    \label{fig:fig-Brain-Intrp0408-Spline-Objs}
\end{figure*}

\begin{figure*}
    \centering
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_04x_Spline_C1.eps}}
    \subfloat[Coil 1]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_08x_Spline_C1.eps}} \\
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_04x_Spline_C2.eps}}
    \subfloat[Coil 2]{\includegraphics[width=3.00in]{C:/TeX/Thesis/Chapter4/Figures/pFig_Brain_08x_Spline_C2.eps}} \\
    \caption{The coil sensitivity maps reconstructed from the experiment of human brain. Cubic spline interpolation was used to coil sensitivity. All images are displayed with $128 \times 128$ resolution. (a) (c) are from $4 \times$ interpolation. (b) (d) are from $8 \times$ interpolation. Most of the artifacts outside of the head are removed. All maps are normalized to $1$. No regularization is used.}
    \label{fig:fig-Brain-Intrp0408-Spline-Coil}
\end{figure*}

\section{Conclusion}
In this chapter, we extended the cubic convolution interpolation, quadratic line search, polynomial approximation of the exponential time function and nonuniform FFT from the single-coil SS-PARSE to the multiple-coil SS-PARSE.

Because the reconstruction complexity is proportional to the number of the receiving coils, the computational speed improved by the fast approach stated in Chapter \ref{chap:2} is more significant in the reconstruction of the multiple-coil system.

The experiments of simulation show that regularization can improve the reconstruction performance. Because we lack the golden standard of human brain, we do not use regularization for the human experiment.

We applied both cubic convolution interpolation and cubic spline interpolation for coils sensitivity. We see similar results for the human experiment. But in simulation, spline can not reach lower error by the quadratic line search, so we have to resort to golden section search that takes much longer time. Even we do not observe this phenomenon in the human experiment, we still believe it is safe to use cubic convolution.  